Turbulent-like Fluctuations in Quasistatic Flow of Granular Media 
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We analyze particle velocity fluctuations in a simulated granular system subjected to homoge- 
neous quasistatic shearing. We show that these fluctuations share the following scaling character- 
istics of fluid turbulence in spite of their different physical origins: 1) Scale-dependent probability 
distribution with non-Guassian broadening at small time scales; 2) Power-law spectrum, reflecting 
long-range correlations and the self-affme nature of the fluctuations; 3) Superdiffusion with respect 
to the mean background flow. 

PACS numbers: 83.80.Fg, 74.80.-g, 45.70.Mg 
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The key role of fluctuations in quasistatic (QS) flow of 
granular media has been noted by several authors refering 
basically to stress fluctuations in time or the inhomoge- 
neous distribution of forces in space Q ^|. Amazingly 
few studies have been reported about the fluctuations of 
particle velocities under homogeneous strain conditions. 
These fluctuations have been observed to occur in a cor- 
related fashion, though their scaling properties have not 
yet been analyzed^, Other recent studies concern 
mainly Couette flows where the strain is localized in the 
vicinity of the inner rotating cylinder]^]. 

Analogies may be drawn with driven multiphase me- 
dia dominated by steric interactions at the microscopic 
scale. Rheologically oriented studies of foams, for exam- 
ple, reveal large-scale cell rearrangements and vortex-like 
structures of the velocity field that control the complex 
flow behavior of foams || However, the (capillary) 
elasticity of cell walls is an important factor that discour- 
ages a closer comparison with systems composed of stiff 
elements such as granular materials. Concentrated sus- 
pensions provide a closer analogy in this respect. There 
exists now convincing evidence based on particle tracking 
that non-Brownian particles in a sheared fluid at small 
Reynolds numbers undergo a (kinematic) diffusive mo- 
tion beyond a concentration-dependent set-in time|8|, 

A much more remote comparison can be made with the 
field of turbulence. There obviously the physics is fon- 
damentally different from that governing granular media. 
Nevertheless, the rich body of work devoted to the statis- 
tical analysis of the fluctuating part of the velocity field 
in turbulence provides a suitable framework that can be 
applied in order to characterize the analogous fluctuat- 
ing part of the velocity field in a granular medium. To be 
more concrete, let us consider the following key aspects of 
turbulent fluctuationsjT^]: 1) Non-Gaussian broadening 
of the probaility density functions (pdf's) of velocity dif- 
ferences as a manifestation of small-scale intermittency, 
2) Multiscale organization of the velocity field reflected 
in its power-law spectrum, and 3) Anomalous diffusion 
pertaining to the Richardson regime. Although deeply 



rooted in fluid dynamics (Navier-Stokes equations, iner- 
tia regime), these scaling features may, in principal, prove 
to be relevant as well within a different physical context 
such as granular flows. 

Following this route, the objective of this paper is to 
show that particle velocity fluctuations in a QS granu- 
lar flow, simulated by means of the molecular dynam- 
ics method, exhibit indeed strikingly similar features. A 
strict analogy with turbulence makes certainly not much 
sense because of a drastically different physics that un- 
derlies these fluctuations, fn particular, inertia effects 
are basically irrelevant in a QS granular flow, and fric- 
tional and hard-core inelastic interactions between parti- 
cles have little in common with molecular interactions in 
a fluid. But, precisely because of these dissimilarities, the 
observed analogy in terms of scaling properties is quite 
nontrivial and it might lead to new insights in both fields. 
In the following, we first describe the simulated granular 
system and our procedures of data analysis. Then, we 
present our main results focussing on the pdf's, correla- 
tions and diffusion, respectively. 

The investigated granular model is a two-dimensional 
assembly of 4000 frictional discs with diameters uni- 
formly distributed between D m i n and D max with 
Dmax — iD m i n . The particles interact through a stiff 
linear repulsive force as a function of mutual overlaps 
and the Coulomb friction law. The equations of motion 
for particle displacements and rotations are integrated by 
means of a predictor-corrector scheme [pd| . 

An accurate evaluation of the statistics of fluctuations 
requires long-time homogeneous and steady shearing. 
However, ordinary wall-type boundary conditions induce 
a pronounced layering effect and the corners enhance the 
local frustrations whereby large strain and stress inhomo- 
geneities arise when the box shape changes. System-size 
inhomogeneities may also occur due to shear localization. 
In order to circumvent such unwanted effects, we used bi- 
periodic boundary conditions following a method similar 
in spirit to that devised by Parrinello and Rahman p^]. 

In our simulations, the gravity was set to zero and 




FIG. 1: A snapshot of particle displacements Ss l with respect 
to the mean background flow. 



a confining pressure was applied along the y direction. 
The width L of the simulation cell was kept constant. 
The displacement field is decomposed into an affine dis- 
placement field 5r % = (Sr x ,Sr y ) and a fluctuating field 
5s l = (5s x ,5s y ) of zero mean {{5s) = 0) carried by the 
particles i. The system is driven by imposing 5r x l — 
Stjry, where 7 is a constant shear rate and 5t is the time 
step. In other words, the k = Fourier mode of the total 
strain is imposed, corresponding to a large scale forc- 
ing. This driving mode was applied on a dense packing 
prepared by isotropic compaction. The height H of the 
packing increases (dilation) in the initial stages of shear- 
ing before a homogeneous steady state is reached where 
5H/H = 5ry '/r y fluctuates around zero. The focus of 
this paper is the field 8s 1 which corresponds to a spa- 
tially periodic motion of the particles with respect to the 
background flow 8r l . 

Although our dynamic simulations involve the physical 
time, the inertial effects are negligibly small and the gran- 
ular texture evolves quasistatically at time scales well 
below 7 -1 . We normalize all times by 7 -1 so that the 
dimensionless time t in what follows will actually repre- 
sent the cumulative shear strain. We will also use the 
mean particle diameter D to scale displacements. As a 
result, the velocities will be scaled by 7D and the power 
spectra in space by (D 2 j) 2 . In our simulations the time 
step was 5t ~ 10~ 7 (in dimensionless units), and more 
than 2.10 7 steps were simulated, corresponding to a total 
strain larger than 2. The solid fraction and the average 
coordination number in the steady state were about 0.8 
and 3.8, respectively. 

It is important to emphasize here that the granular 
nature of our system does not allow to apply exactly the 
same procedures of data analysis as in turbulence. Tur- 
bulence studies focus mainly on velocity differences 5v 
measured at a fixed point of a fluid over a time interval r 
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FIG. 2: The pdf 's of the y-components of fluctuating veloc- 
ities for two different integration times: 10 -3 (broad curve) 
and 10 _1 (narrow curve). The latter is fitted by a Gaussian. 
The error bars are too small to be shown. 

or between two points separated by a distance r. In con- 
trast, the particle-scale granular motion involves a dis- 
crete displacement field that is carried by individual par- 
ticles. Thus, our natural framework is Lagrangian rather 
than Eulerian. These differences are certainly important 
for a strict one-to-one comparison, but here we basically 
consider turbulence as a reference field from which we 
extract tools to characterize granular fluctuations. 

Another distinctive feature of granular flow is that, 
due to collisions, the velocities are discontinuous in time. 
As the positions are better behaved, we characterize the 
fluctuating motions of the particles by "tracer" velocities 
defined from particle displacements Ss l by 

1 r t+T . 

v%t + r) = -J 5s\t') dt' (1) 

where r is the time resolution. Since we are concerned 
with steady flow, the statistical properties of v (pdf's, 
correlations) are independent of t (though, as shown be- 
low, they crucially depend on r). Hence, accurate statis- 
tics can be obtained by cumulating the data from differ- 
ent time slices of a single simulation running for a long 
time. 

Fig. [l] shows a snapshot of fluctuating velocities v l 
for a short time lag r = 10~ 7 . We see that large-scale 
well-organized displacements coexist with a strongly in- 
homogeneous distribution of amplitudes and directions 
on different scales. Eddy-like structures (though without 
the singular vorticity concentrated in the core of these 
eddies) appear quite frequently, but they survive typi- 
cally for strains r less than 10~ 3 . After such short times, 
large-scale eddies break down and new (statistically un- 
corrected) structures appear. This behaviour is radically 
different from turbulence where the eddies survive long 
enough to undergo a significant distortion due to fluid 
motion |l0[]. 

We consider now the pdf's of the components v x (t, t + 
t) and v y (t, t + r) as a function of time resolution r. The 
pdf's of v y are shown in Fig.^J for a short integration time 
t = 1CP 3 , and for a long integration time r = 10 _1 . The 
pdf has changed from a nearly Gaussian shape at large r 
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FIG. 3: Flatness F of the distribution of velocity fluctuations 
as a function of the integration time r. 



FIG. 4: Averaged power spectrum of the x and y components 
of the fluctuating velocity field over one time step for one- 
dimensional cross sections along the mean flow. 



to a non-Gaussian shape with broad stretched exponen- 
tial tails extending nearly to the center of the distribution 
at small r. We found no simple form allowing to fit the 
non-Gaussian pdf over the whole velocity range. 

In order to characterize this non-Gaussian broadening 
of the pdf 's as a function of r, we calculated the flatness 
F = (Vy)/(Vy) — 3, which is zero for a Gaussian distri- 
bution and 3 for a purely exponential distribution. The 
values of F as a function of r, shown in Fig|j|, are con- 
sistent with zero at large r (r > 0.2) and rises to 5 for 
our finest time resolution (r = 10~ 7 ). A strictly similar 
behaviour was observed for the component v x . 

The broadening of the exponential tails of the pdf 's at 
increasingly smaller scales is a hallmark of fully devel- 
opped turbulence (for velocity differences) fill . It is at- 
tributed to the phenomenon of intermittency, i.e. strong 
localized energy transfers at small scales. In a QS gran- 
ular flow, the basic physical mechanism underlying the 
fluctuations is the mismatch of the uniform strain field 
applied at the boundaries or in the bulk, with mutual 
exclusions of the particles. As a result, the local strains 
deviate from the mean (global) strain. The observation 
of a transition toward a Gaussian distribution for large 
time lags is a sign of loss of correlation and/or exhaus- 
tion of large fluctuations in the increment of displace- 
ment which occur at different times. Unfortunately, the 
rich multifractal scaling of velocity fluctuations is out of 
reach within the present investigation due to demand- 
ing statistics jl^]. The analogy with turbulence, however, 
suggests further study along such routes. 

In order to quantify the extent of these correlations, 
we estimated the power spectrum E of velocity fluctu- 
ations both along and perpendicular to the flow and at 
different times. The Fourier transform was performed 
over the fluctuating velocity field defined on a fine grid 
by interpolating the velocities from particle centers. The 
power spectra were quite similar along and perpendicular 
to the flow, and for different snapshots of the flow. The 
averaged one-dimensional spectrum is shown in Figji| It 
has a clear power-law shape k~^ ranging from the small- 
est wavenumber k = D/L, corresponding to scale L, up 
to a cut-off around k = 0.5, corresponding to nearly two 
particle diameters. The exponent is (3 ~ 1.24 ~ 5/4 over 



one decade (to be compared with the exponent 5/3 as 
a hallmark of 3D turbulence for the spectrum of veloc- 
ity differences). This means in practice that the fluctu- 
ating velocity field is self-affine with a Hurst exponent 
#=03-l)/2 = O.12||. 

Due to the peculiar behavior of the velocity field (dis- 
continuous in time), one might expect that the power 
spectrum is sensitive to the time resolution r. However, 
we checked that the value of /3 is independent of r. It is 
also noteworthy that, the presence of long range correla- 
tions in displacements, reflected in the value of f3, is in 
strong contrast with the observed correlation lengths of 
nearly 10D for contact forces^]. 
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FIG. 5: Diffusive trajectory of one particle with respect to 
the mean background flow for a cumulative shear strain of 2. 

The long-time behaviour may be studied by consider- 
ing the effective diffusion of the particles. Normal dif- 
fusion implies that the root-mean-square (rms) relative 
displacements A in a given direction varies in propor- 
tion to the square root of time. In 3D fluid turbulence 
at high Reynolds numbers, the long-time pair diffusiv- 
ity of suspended particles is anomalous, following the 
Richardson law A cx r 3 / 2 . This superdiffusion law reflects 
the Kolmogorov-Obukhov velocity spectrum v\ cx A 1 / 3 , 
where v\ — A(r)/r is the characteristic velocity differ- 
ence over a distance A|l(| |lfj| . 

We analyzed the kinematic diffusion of single particles 
in our QS granular flow. One example of a single particle 
trajectory with respect to the background strain is shown 
in Fig.^. We see that the fluctuating displacement is of 
the order of the mean particle diameter for a strain of the 
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FIG. 6: Root-mean-square relative displacements A along 
x and y directions as a function of time r fitted by a power 
law of exponent 0.9 (straight solid line). As a guide to the 
eyes, a power law of exponent 1/2, corresponding to normal 
diffusion, and the line A oc r are shown as well (dashed lines). 
The error bars are too small to be shown. 

order of unity. Fig.|^ shows the rms relative displacement 
A(t) of all particle pairs initially in contact, as a function 
of time along x and y directions. This clearly corresponds 
to a superdiffusion behaviour, 

A oc T a , (2) 

with a ~ 0.9 for both components over nearly 3 decades 
of strain. Particle self-diffusivities exhibit a similar law. 
Since large-scale structures are short lived, anomalous 
diffusion scaling in our granular system can not be solely 
attributed to velocity correlations. This behaviour re- 
veals, above all, the long-time configurational memory of 
a granular medium in QS flow. 

We note that Fig.|] shows no anisotropy for the diffu- 
sion. In fact, the steric exclusion effects dominate over 
the large-scale strain field for small diffusive displace- 
ments (one particle diameter, corresponding to 2% of the 
cell size in our simulations), and hence the anisotropy 
may be weak at such scales. Moreover, we observe no 
crossover to normal scaling within the investigated strain 
range. We cannot exclude that for larger strains, when 
two particles are more widely separated, a normal diffu- 
sion law is recovered. However, this would require ex- 
tremely long computations, and, experimentally, strains 
of the order of unity are already above standard tests. 

In summary, we analyzed fluctuating particle displace- 
ments with respect to the background quasistatic shear 
flow in a model granular medium. These fluctuations 



were shown to have the following scaling charactersitics: 
1) The pdf's undergo a transition from streched expo- 
nential to gaussian as the time lag is increased; 2) The 
power spectrum (in space) of the velocity field obeys a 
power law, reflecting long range correlations and the self- 
affine nature of the fluctuations; 3) The fluctuating dis- 
placements have a superdiffusive character. These ob- 
servations contradict somehow the conventional wisdom 
which disregards kinematic fluctuations in macroscopic 
modelling of plastic flow in granular media. Several ba- 
sic aspects of quasistatic granular flow (elementary rep- 
resentative volumes, mean field approximation, memory 
effects, mixing) are thus concerned by these findings |l7[ . 

There appears an evident analogy with the scaling fea- 
tures of turbulence that was also discussed throughout 
this paper. Of course, this analogy does not imply that 
a QS granular flow can be considered as turbulent in the 
standard sense of fluid dynamics. In particular, because 
of the fundamentally different origin of granular fluctua- 
tions, a direct reference to the underlying physics of tur- 
bulence can be misleading. But, the observed analogy 
in terms of scaling characteristics is consistent enough 
to upgrade kinematic fluctuations in quasistatic granular 
flows to the rank of a systematic phenomenology which 
could be coined by the term "granulence" as compared 
to "turbulence" in fluids. 

Interestingly, this analogy works with three- 
dimensional turbulence although our data concern 
a two-dimensional granular flow. The energy cascade in 
turbulence is governed by inertia, and two- and three- 
dimensional systems do differ significantly in this respect. 
In quasistatic granular flow, the fluctuating velocity 
field is a consequence of the geometrical compatibility 
of the strain with particle arrangements, and dissipation 
is mainly governed by friction at the particle scale. The 
difference between two- and three-dimensional systems 
may thus be less crucial, but it was not investigated in 
the present work. Quite independently of its physical 
origins, this analogy is suggestive enough to be used 
as a promising strategy towards a refined probabilistic 
description of granular flow in the plastic regime. 
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